namespace hms::swe {

template<typename MeshType>
template<BlockLocation loc>
void SolverThread<MeshType>::computeFriction(
	const Block<loc>& block,
	Array3XsMap& fluxMap,
	const Ref<const ArrayXXs>& dBlock,
	const Ref<const ArrayXXs>& vBlock
){
	scalar dt { solver().time().step() };
	/* allocate memory for helper variables */
	// Array1Xs vNorm, chezyFac,fricFac;
	// vNorm.resize(1, block.nx() );
	// chezyFac.resizeLike(vNorm);
	// fricFac .resizeLike(vNorm);
	/* don't allocate, but map memory (30% faster for overall friction computation) */
	assert( m_fricBuf.cols() >= 3 && m_fricBuf.rows() >= block.nx() );
	Array1XsMap
		vNorm    { m_fricBuf.col(0).data(), block.nx() },
		chezyFac { m_fricBuf.col(1).data(), block.nx() },
		fricFac  { m_fricBuf.col(2).data(), block.nx() };
	assert( chezyFac.data() >= vNorm.data() + vNorm.size() );
	assert( fricFac.data() >= chezyFac.data() + chezyFac.size() );

	/* Friction source term with preconditioning */
	/* See Simons (2019), Eq. 3.109. */
	std::visit( [&](const auto& nField){

	block.forEachRow( [&](Index y){

		auto fluxes { block.rowMap_mem(fluxMap, y).template bottomRows<2>() };
		auto d { block.template rowMap<1>(dBlock, y) };
		auto v { block.template rowMap<2>(vBlock, y) };
		auto n = [&](){
			if constexpr ( remove_qualifiers<decltype(nField)>::isUniform ){
				return
					block.template writeBlock<1>( nField.values() )
					.col(y).transpose();
			} else {
				return block.template rowMap<1>(
					block.template writeBlock<1>( nField.values() ), y
				);
			}
		}();

		vNorm = v.matrix().colwise().norm().array();
		/* with
		 *         C = d ^ 1/6 / n
		 * <=>   C^2 = d ^ 1/3 / n^2
		 * <=> 1/C^2 = n^2 / d^1/3
		 * */
		chezyFac = dt * phys::g * n.pow(2) / d.pow(1./3);
		fricFac = chezyFac * vNorm;
		chezyFac /= d;
		
		auto vNormRepl { vNorm.template replicate<2,1>() };
		auto fricRepl  { fricFac.template replicate<2,1>() };
		auto chezyRepl { chezyFac.template replicate<2,1>() };

		fluxes = ( d < SolverType::minWaterDepth || vNorm == 0 )
			.template replicate<2,1>().select(
				fluxes,
				( fluxes - fricRepl * v ) /
				( 1 + chezyRepl * ( vNormRepl + v.pow(2) / vNormRepl ) )
			);
		
		/* the equation boils down to
		 *                 g
		 *     f_0 - dt * --- * |v| * v
		 *                C^2
		 * f = ----------------------------
		 *                 g              v^2
		 *     1 + dt * ------- * ( |v| + --- )
		 *              C^2 * d           |v|
		 * 
		 * with f_0 being the momentum fluxes from the flux scheme
		 * and f being the adjusted fluxes.
		 * */
	});

	}, solver().n() );
}

} // namespace hms::swe